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COMPUTATIONAL  ASPECTS  OF  THE  PROBLEM  OF 
OPTIMAL  FLIGHT  AS  A  BOUNDARY  PROBLEM 

V.  K.  Isayev  and  V.  V.  Sonin 
(Moscow) 

I.  Work  [l]  proposed  a  modification  of  the  Newton  method  for 
solving  boundary  problems  for  ordinary  differential  equations.  In 
this  article  we  examine  one  example  of  applying  this  method  to  the 
problem  of  optimal  motion  of  a  point  of  variable  mass  M(t)  In  a 
central  gravitational  field  with  a  limited  expenditure  of  power. 

For  convenience  let  us  briefly  reproduce  the  derivation  of  the 
basic  relationships  (for  more  detail,  see  [2j).  Let  the  allowable 
velocity  £  of  the  discharged  mass  be  limited  0  <  S  c^  <  00 ,  and 
let  the  power  of  the  reactive  jet  N  =  -  Me  /£?  lie  within  the  limits 
0  s  N  ^  N_„  for  all  values  of  c.  The  flat  motion  of  a  mass  point 
in  vacuum  in  a  Newtonian  gravitational  field  is  described  by  the  system 


u  =» 

Nut  cos  9 

iix 

me 

(** +  »*)*/•  ’ 

Nut  sin  <p 

Mf 

V  tiT- 

me 

(*,  +  y,),/* 

Kttt 

m  a*  — r*, 

c2 

where  u  j  -  N/N^,  N  =  ?N1aw/M(  0 ) ,  m(  t )  =  M(t)/M(0);  x,  y  -  the  Cartesian 
coordinates  of  the  point  of  variable  mass,  u,  v  -  components  of  Its 
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velocity,  <p  —  the  angle  between  the  Ox  axis  and  the  direction  of  the 
thrust  vector.  We  shall  also  search  for  the  programs  of  ~hange  in 
power  u 1 ( t )  (0  s  u  (t)  =  l),  the  velocity  of  the  reactive  jet  c(t) 
and  the  direction  of  the  thrust  vector  <p(t),  which  during  the  time 
t  -•  T  move  the  material  point  with  an  initial  ma3s  M(0)  from  the 
position  x (  0 )  ---  x°,  yi'O)  -  u°,  u(0)  =  v°  to  the  position  x(fr)  =  x1 , 
y(T)  =  y1 ,  u(T)  =  u1 ,  v(T)  =  v1  under  conditions  of  minimal  expenditure 
of  mass  or  equivalent  conditions  of  maximal  final  mass. 

2.  To  solve  the  stated  problem  with  the  aid  of  the  maximum 
principle  of  L.  S.  Pontryagin,  let  us  compose  the  function 


11  —  pu 


[ 


Nui  cos  <p 
me 


_ f*f _ 1 

,  Nu\ 

+  +  —  Pm~~, 

c * 


+  p  l> 


Nii  sin  9p 
me 


py  "I 


where  pt(i  ~  u,  v,  x,  y,  m)  are  the  conjugate  variables  that  satisfy 
the  system  of  differential  equations 

Oil 

P(==— —  (i  =  u,v,x,y,m).  (2) 

The  condition  of.  the  maximum  of  the  function  It  relative  to  hie  control 
functions  (u^.r,  cp)  for  the  given  bounds  gi’-s  values  of  the  max- 
-  optimal  equations,  that  are  introduced  in  Table  1. 


Table  1 


Valuta  if  the'Vftlaal  ofeiatMM 
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m  cmin 

l 

emln 

V?u  +  c  2pm  «m.« 

!  t 

t*.  _ 

V  p*«  +  p*. 

m 

3 

-Pm<‘  2^r, 

i 

em*t 

o 

i. 


FTD-HT-67-44  7 


Putting  the  values  of  the  optimal  equations  found  from  the  table 
into  the  initial  (l)  and  conjugate  (2)  systems,  we  reduce  the  vari¬ 
ational  problem  formulated  above  to  a  boundary  problem  of  the  10-th 
order  with  these  conditions: 

“(0)“u®  v{0)  =r  v°,  X(0)=x“,  i/(0)  =  y°;  m*-i;  (3) 

u(T)  =  «i,  v(T)  =  t,*,  X(T)  =  x«,  y(T )  =  y\ 

Pm(T)  —  prr}  —  — 1.  (  4  ) 

The  variational  problem  was  reduced  to  calculating  p  (0),  p  (0), 

PjtiO.  p,(0),  p.(0)  from  conditions  (4),  which  were  fulfilled  on  the 
right  end. 

We  shall  consider  that  the  values  u(T),  v(T),  x(T),  y(T),  p  (T) 
are  well  approximated  by  the  continuous  and  differentiable  functions 
of  the  nought  initial  conditions.  Let  there  be  an  approximation 
P^O)  =  Pt°  (1  =  u,  v,  x,  y,  m)  that  corresponds  to  the  solution  of 
the  system  of  equations  (l),  (?)  that  satisfies  conditions  (3),  but 
does  not  satisfy  conditions  (4).  If  there  are  small  perturbations 
of  the  sought  approximation 

/',(())  +  W  (<*«.*.*,„,»),  (b) 

because  of  which  the  Cauchy  problem  for  the  system  of  equations  (l), 

)  with  the  initial  conditions  (3),  (5)  satisfies  conditions  (4), 


“ ( T'  pu° +  P"*  +  W,  p*°  4-  ftp*0,  p¥*  +  6 p„\  Pm*  -f  tPm •)  „  u.t 

v(T,  pu*  +  ftp.0,  pv°  +  ftp,,*,  Pt»  +  jPioi  pyo  +  ^  pm,  +  6pme)  v,f 
*{T,  P«°  +  ftp.0.  p,°  +  ftp.*,  p*  4*  ftp/,  pS  +  tpj,  pm*  +  tPmo)  „  x,f 
V(T,  P*°  +  ftp.9,  p. 0  4-  ftp.*,  p,*  4-  ftp/,  p/  4-  «P*#,  Pm*  4-  ftp*0)  -  p‘, 

Mr,  p„*  4-  ftp-*,  p.*  4-  ftp,*,  p,*  4-  ftp,*,  P»*  4-  ftp/.  Pm*  4-  fipm*)  -  Pm*, 

then,  by  disregarding  the  terms  higher  than  the  first  order,  we 
can  write  (6)  in  the  form 
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du  du  du  du  du 

_4p..+_4,..+_s?,.+_w+__w. 

—  M1  —  U(T,  Pu\  Pu®,  p*°:  Py°,  Pm°) , 
dv  dv  dv  du  dv 

dp?bPu9  +  +  T3»*  +  XT. 6*#  +  TrU>bPm  " 


0p*®^'  '  dPy°  r"  '  Qpm*  r 
•  i>*  —  y(f,  pu#,  pe#,  p,*,  P„#,  pm®), 


dx  dx  dx  dx  dx 

- 6Pu*  4 - flp*°  4- - ftp*®  + - .ftp/  4 - „6p** 

0p„®  P  T  dp.®°P*  T  dP/P*  T  <>Pv  ™  dPm°  .  . 

■»**--*( r,  pu®,  p,#,  p«®,  Pv®,  Pm*) ,  '  '  ' 

=3‘Vl  —  V(T,  Pu®,  Pu®,  Px®,  P»®.  Pm®)  , 

^Pm «  ,  ,  ^Pm  ,  4  ,  3pm  #  dpm  flpn»  ,  a _ 

«^4p-  +«w+  I*?4'’'  +  4p‘  +  'ST? 4p-  " 

“Pm1-  Pm  ( T,  pu°,  Pu®,  p,®  ,P„®,  Pm®) . 

In  order  to  use  system  (7),  we  must  calculate  the  values 

U(T,  Pu®,  P„®,  Pi®,  P»®,Pm°) ,  W(r,  Pu®,  P»®,  Pi®,  Pi,®,  Pm®),  *(7\  Pm®,  P»®.  Pi®,  P»°-  P«°).  P(7’>  P“®<  Pv#»  P*#,  P»#> 

Pm®),  Pm(7’,pu®,pn#,p,#, p/,pm°),  that  stand  in  the  right  parts  of  (7)  and  we 
must  calculate  the  matrix  of  the  derivatives. 

We  shall  obtain  the  first  by  solving  the  Cauchy  problem  for  the 
system  (l),  (?)  with  the  initial  conditions  (d)  and 

Pu(0)  =  Pu®,  P«(0)  ~  Pu®,  Px(0)  —  Pi®,  P»(^)  ~  P»®> 

Pm{0)  =  pm®.  (  3  ) 

To  calculate  each  column  of  the  matrix  of  derivatives,  we  must  solve 
the  system  of  equations  in  variations  or  (additionally)  we  must  solve 
the  Cauchy  problem  for  sysi.em  (l),  (?). 

For  example,  the  first  column  of  the  matrix  Is  determined  from 


the  formulas 


du  _  U(T>  P“*  4-  Apu®,  Pu*.  Pi®,  p»®,  Pm*)  -  U(T,  Pu®,  p.®,  p.®,  P»®,  Pm,®) 


*Pm  —  Pmtf,  Pu*4-Ap,*.  P,®,  pi®,  p/.  Pm®)  ~Pm(T,  PJ,  p,\  p««,  p,»,  pw«)  (  V  ) 

dpu*  Ap.# 

In  order  to  r.se  formulas  (y),  we  must  solve  the  Cauchy  problem 
for  the  system  (1),  (?)  with  the  initial  given  conditions  (j)  and 
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P.(0)-p.'  +  ip.-.p.(0)=Pl.,  P,(0)=p,.,  p,(0)=p,..  ,„«,)=  p„. 

The  remaining  4  columns  of  the  matrix  can  be  calculated  in  a 
similar  manner. 

To  solve  the  Cauchy  problem,  we  used  the  Runge-Kutta  method 
with  automatic  selection  of  the  step  with  respect  to  accuracy. 
During  calculation  of  the  derivatives  we  took 


\Pt°\>  1; 

| /’.“!<  1  (i  =  u,v  x,i/,m). 


After  ier.ermining  the  numerical  values  of  the  parameters  of  system 
v  7 ) .  Let  us  find  the  increments  dpt°(i  -  u,  v,  x,  y,  m)  from  it. 

Taking  pt°  +  6pl°  (i  =  u,  v,  x,  y,  m)  as  the  initial  point,  let  us 
repeat  the  procedure  of  determining  the  parameters  of  the  "new"  system 
.7)  and  calculating  the  "new"  dp^0  . 

However,  the  clasical  Newtonian  method  introduced  above,  as  it 
applies  to  the  examined  problem,  usually  does  not  converge.  To  Improve 
the  convergence,  according  to  the  terminology  of  work  [l],  let  us 
introduce  the  "miss  function,"  i.e. ,  the  magnitude  that  characterises 
the  quality  of  fulfillment  of  conditions  ^4).  (In  the  given  example 
the  distance  between  the  actual  end  of  the  phase  trajectory  and  the 
point  (4)  was  taken  as  this  magnitude.  The  use  of  this  function 
allows  us  to  correct  the  found  increments  dp^0  (1  -  u,  v,  x,  y,  m), 
which  improves  the  convergence.  ) 

For  this  purpose,  after  determining  dp^  (i  =  u,  v,  x,  y,  m)  let 
us  calculate  the  functions 

u(7\p.  +  alps',  p,®  -f  alp,',  p,'  -f  alp,*,  p,*  +  o bp,*,  pm*  +  o lpm*)  -* 

-■(r,p*+«ep®).  (10) 

. . .  * . t  ..•  # 

Pm(T,  p,'  +  alp,*,  p,'  +  alp,',  p,*  +  alp,*,  p,*  -f  alp,*,  pm*  +  alpm*)  — 

“  Pm(T ,  p®  +  O ftp*), 


solving  for  the  series  of  fixed  values  a  the  Cauchy  problem  for  systems 
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(1),  (2)  with  the  initial  conditions  (3^  and  p„  (0)  =  p*#  -f  adp„#,  p.(0)=»pe44- 
+  a flpe\  p,(0)  =  p,#  +  a6p«°,  Py (0)  =  p„°  +  a6pv#,  Pm{0)  *=  pm*  +  aty*4.  From  the  values 

(10)  let  us  find  the  "miss  function"  <t>(a) 

©i(a)  =  [u(r,  p'  +  abp0)  -u*]l+  [viT.pP  +  aty)  -  «>'J* + 

4-  [x{T,  p°  +  a6p°)  -i*]*+  [p(r,  P®  +  a6p#)  -  V*)*  + 

4- [p-n(f,  p°  4- a6p°)  —  Pm1]1 

and  the  value  a*  that  corresponds  to  the  minimum  $>(a)  in  the  interval 

A 

0  <  a  §  1.  From  the  preliminarily  calculated  values  of  4>(i)  and 
let  us  find  the  parabola  that  approximates  0(a): 

©(a) «  A(a),-  ©(0)-  [  30(0)—  4©  (  A)  4-  ©(i)  ]  a  4- 

4*  2  [  ©  (0)  —  2©  ^  A  )  4-  ©  ( i )  j  a*. 

The  parabola  Lj(a)  has  one  extremum  at  the  point 

3©(0)-4©(A)+©(1) 

Of  WLM  — _ _ _ _ _ _ 

4[©(0)-2©(Aj+©(i)J. 

Let  us  examine  the  following  separately: 

1)  the  case  of  the  minimum  ©(0) -2©(  ^)  4-©(l)  >0:  •)  0,^1, 

b)  ^  <  a*  <  i,c)  a,  <  0; 

2)  the  ca3e  of  the  maximum  ©(0)  —  2©{,i )  4-  ©(I)  <  0:  a)  a,  >  \,  t)  0  <  n,  <  l, 

O)  a,  ^  o. 

In  cases  l,n)  and  , c )  and  with  the  additional  condition  that 
1 )  <  t>[0)  Is  fulfilled,  also  in  case  T,b))  we  should  set  a*  =  1; 
in  case  l,b),  we  should  take  a*  -  a#;  In  remaining  cases  we  should 

i 

repeat  the  process  in  the  lesser  interval  0  %  a  §  i.e.  ,  we  snould 

4  4  4 

calculate  fr~m  the  magnitudes  of  i(0),  and  l\-£)*we  should 

construct  the  parabola  L^a)  that  approximates  t\a)  In  the  interval 
0  i  a  S  x»  etc.  If  in  this  case  we  do  not  manage  to  find  the  value 
a*,  we  should  calculate  ©(. g-) ,  construct  tne  polyr.ominal  L»(a)  for  the 
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interval  0  ^  a  %  etc. ,  until  the  process  leads  to  a  certain  value 
of  a*.  After  a*  is  found,  the  iterative  process  is  renewed  from 
point  ( 3 )  and 

MO)  =  M  +  afipj,  pv( 0)  =  pt°  |  «*Apt.® 

Px(0)  =  px°  +  a*6px#,  p„(0)  =  p,#  +  a*fip/, 

Pm(0)  =  pm°  -(-  a*fipm». 

5.  Let  us  make  a  number  of  comments  connected  with  the  realization 


of  the  described  algorithm  on  a  computer. 

With,  this  method  a  Cauchy  problem  is  solved  for  a  system  of  10-th 
order  ejections  v  1 ) ,  ( 7 )  for  lifferent  initial  conditions  at  each 
step  of  the  iteration  not  less  than  eight  time-  , six  times  to  obtain 
the  right  parts  and  the  matrix  of  rhe  derivatives  of  system  (7)  and 
twice  to  d  etc ’--mine  1\k)  and  11)),  a  syst  em  of  linear  5-th  order 
equations  Is  solved,  a  series  of  logical  operations  is  carried  out. 


As  •omput.ational  practl  -e  shows,  at  least  y0#  of  the  machine 
time  is  spent  on  solving  the  Cauchy  problem.  Therefore,  to  reduce 
the  machine  time  it  Is  expedient  at  the  start  of  the  iterative  search. 


whet;  the  "miss"  with  respect  to  the  boundary  -on Jit  ions  1?  still 
comparatively  great,  to  solve  the  Cauchy  problem  with  a  relatively 
largo  constant  step  or  with  a  relatively  large  error  in  the  step 
,  during  calculation  with  automatic  selection,  of  the  step),  decreasing 
the  st  ep  at:  1. 1)  is  reduced  ;  or,  correspondingly,  the  error  that  is 
all  ok  die  in  the  step). 

While  the  accuracy  of  the  calculation  will  not  Increase  in  the 
future,  the  computational  error  will  become  comparable  wltn  the 
magnitude  of  the  miss  and  the  process  will  sense  tc  converge.  Hence 
It  is  evident  that  it  13  necessary  to  control  the  accuracy  of  the 
solution  to  a  system  of  differential  eouatlone  in  the  Interval  [0,  T). 


*:■<*•*«*•. 


We  used  the  following  method  :  system  (l)  and  .2)  has  two  first 
integrals  (see  [?])  H  -  const  and  M  =  pnv  -  pvu  +  pxy  -  pyx  =  const. 

On  the  ends  of  the  integration  interval  we  calculated  the  values 
M( 0 ) ,  H(T),  M(T),  after  which  the  larger  of  the  magnitudes 

\H(0)-~II(T)\  jJ/(0)-A/<7’)| 

|//(0)|  •  i  a/to,  r 

took  into  account  the  relative  error  e  of  the  solution  in  the  interval 
[ 0,  T].  If  the  relative  error  became  less  than  e,  then  the  accuracy 
of  calculation  on  the  following  iteration  step  increased. 

4.  Let  us  examine  a  specific  example  of  the  solution  of  the 
boundary  problem  for  the  following  values  of  the  parameters  and  the 


if* 


boundary  conditions  : 


N  17.506613, 

cmln  —  0.6, 

Cipax  — —  lUGO, 
li  =  39.4784176, 
T  «  0.9 


„o  _ 


0.00, 


0.28318531, 
=  i.00, 
y°  —  0.00, 

7T.°r=  1.00, 


=  0.00 

c*  =5.090159^1, 
i;  =  1.52369100, 
y'  =  0.00, 

Pm1  =  -1.00. 


As  the  initial  approximation  of  p°  we  selected  the  point  pu#  =  0.119420502, 

P„°  =  0.182035475;  p*°  =  1.16399559,  Pv°  =  0.63037049,  pm“  =  0.504237974. 

The  course  of  iteration  is  easily  followed  from  Table  2  and  Pigs. 
1—3  (given  are  the  values  of  the  phase  coordinates,  the  variables 
that  are  connected  to  them,  the  miss  function  0  and  the  first 
integrals  V  and  M  when  t  =  0  and  t  =  T  depending  upon  the  number  of 
iterations  k). 

Figure  4  introduces  certain  coefficients  of  the  matrix  of  the 
system  of  equations  (5)  as  a  function  of  the  parameter  k  (for  conveniene 


* 

The  use  of  the  first  Integrals  to  evaluate  the  accuracy  of  the 
solution  of  a  system  of  ordinary  differential  equations  on  the  interval 
was  used  previously  by  V.  A.  Yegorov  (see,  for  example,  [3]). 

The  problem  of  flight  between  a  circular  orbit  of  earth  and  Mars 
that  is  optimum  with  respect  to  expenditure  of  the  working  mass  (the 
heliocentric  angle  of  flight  f  -  3600,  duration  T  =  0.9  year)  reduces 
to  the  given  boundary  problem. 
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discrete  values  on  Figs.  1  —  4  are  connected  fcy  a  smooth  curve). 

The  intervals  between  whole  neighboring  values  of  the  parameter  k  on 
Fig.  3  conditionally  correspond  to  the  full  lengths  of  the  step  in 
the  classical  Newtonian  method,  and  the  internal  points  describe  the 
behavior  of  the  miss  function  in  the  process  of  "step  subdivision. " 

The  dotted  line  represents  the  interpolation  Lagrangian  polynominal 
L  (a)  (in  this  variant  it  is  not  necessary  to  calculate  L2(a)).  From 
Figs.  1—4  it  is  evident  that  there  is  a  divergence  of  the  classical 
Newtonian  method  in  the  second-third  steps  of  the  iteration  (the 
change  in  the  character  of  the  behavior  of  the  matrix  elements  shown 
or.  Fig.  4  3s  evidently  connected  with  this),  and  when  k  >  7  the  values 
of  all  the  magnitudes  stabilize.  The  materials  in  Table  2  reflect  the 
role  of  the  control  of  calculation  accuracy  in  the  process  of  solving 
the  boundary  problem. 

In  the  11-th,  12-th,  13-th  and  ±4-th  steps  of  iteration,  when 
the  error  In  calculating  the  first  integrals  becomes  comparable  to 
T(k,  a*),  the  solution  of  system  (i),  (2)  is  repeated  for  the  same 
initial  data,  but  with  a  greater  accuracy  in  the  step,  A  comparison 
of  the  first  integrals  shows  that  the  required  degree  of  accuracy 
s  =  10” 7  Is  retained  on  the  resultant  trajectory. 


Submitted 
24  February  1964 
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